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Abstract 

We present a new approach to the problem of alternating signs for fermionic 
many body Monte Carlo simulations. We demonstrate that the exchange of 
identical fermions is typically short-ranged even when the underlying physics 
is dominated by long distance correlations. We show that the exchange 



process has a maximum characteristic range of \/2(l — f)(3h lattice sites, 

where (3 is the inverse temperature, h is the hopping parameter, and / is 

the filling fraction. We introduce the notion of permutation zones, special 

regions of the lattice where identical fermions may interchange and outside of 

which they may not. Using successively larger permutation zones, one can 

extrapolate to obtain thermodynamic observables in regimes where direct 

simulation is impossible. 
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I. INTRODUCTION 



One of the most important challenges for quantum field theory and many body simu- 
lations is the study of light fermion dynamics. Today most fermion algorithms used in 
lattice field theory are based on pseudofermion methods, which calculate the contribution 
of fermions indirectly by means of an effective non-local bosonic action. The most popular 
approach, Hybrid Monte Carlo, deals directly with non-local actions [|I| . Molecular dynam- 
ics is used to propose new configurations and a Metropolis criterion is used to accept or 
reject updates. Although the advantages of pseudofermion methods are substantial, there 
are still important motivations for doing simulations with explicit fermions. The primary 
concern is the question of what fermions actually do in lattice simulations. In lattice gauge 
theory, for example, one might be interested in seeing the effects of confinement on quark- 
antiquark pair separation or the relation between chiral symmetry breaking and light quark 
mobility. These questions are most easily addressed in a local field framework which keeps 
the fermionic degrees of freedom explicit. 

In order to do get anywhere with explicit fermion simulations one must of course address 
the fermion sign problem. In this paper we introduce a new approach to the sign problem 
which has applications to quantum simulations at finite temperature. Unlike the fixed-node 
approach our method makes no assumption about the nodal structure of the eigen- 

functions. It also is not a resummation technique, the underlying principle powering the 
meron-cluster algorithm |J and diagonalization/Monte Carlo methods ||[7|. The approach 
we introduce here is based on the observation that in most finite temperature simulations 
fermion permutations are short ranged. This holds true even for systems with massless 
modes and long distance correlations, as we demonstrate with two examples. In this paper 
we focus on the application of the zone method to simulations with explicit fermions. These 
methods have recently been applied to study chiral symmetry breaking in massless quantum 
electrodynamics in 2 + 1 dimensions ||. The extension to pseudofermion methods and, 
in particular, applications to Euclidean lattice gauge theory will be discussed in a future 
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publication. 



II. WORLDLINES 

We begin with a brief review of the worldline formalism ||. We introduce the basic 
ideas in one spatial dimension before moving on to higher dimensions. Let us consider a 
system with one species of fermion on a periodic chain with L sites, where L is even. Aside 
from an additive constant, the general Hamiltonian can be written as 



a i a i+l | + 5Z °i a i a i- 



H = -hJ2 [4+1°* + 

i i 

Following H we break the Hamiltonian into two parts, H e and H , 



H, 



i even/odd 



We note that H = H c + H . 

We are interested in calculating thermal averages, 



(O) 
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where (3 = (ksT) -1 . For large N, we can write 
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where 
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S e/o = exp(-^if e /o)- 



Inserting a complete set of states at each step, we can write Tr [exp(—/3H)] as 



(z \S o \z 2 n-i) ■■■ (zi\S c \z ) . 



(5) 
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In Fig. 1 a typical set of states \zq) , |^2iv-i) are shown which contribute to the sum 
in (ffl). We call such a contribution a worldline configuration. The shaded plaquettes 



represent locations where S c or S D acts on the corresponding local fermionic state. The 
classical trajectory of each of the fermions can be traced from Euclidean time t — to time 
t = 13. In the case when two identical fermions enter the same shaded plaquette, we adopt 
the convention that the worldlines run parallel and do not cross. With this convention the 
fermion sign associated with Fermi statistics is easy to compute. The worldlines from t = 
to t = P define a permutation of identical fermions. Even permutations carry a fermion 
sign of +1 while odd permutations carry sign —1. The generalization to higher dimensions 
is straightforward. In two dimensions, for example, exp(— j3H) takes the form 

N 

(7) 



exp 



« (sys*szs*) N . 

The sum over all worldline configurations can be calculated with the help of the loop 



algorithm | lQfl . At each occupied/unoccupied site, we place an upward/downward pointing 
arrow as shown in Fig. 2. Due to fermion number conservation, the number of arrows 
pointing into a plaquette equals the number of arrows pointing out of the plaquette. New 
Monte Carlo updates of the worldlines are produced by flipping the arrows which form closed 
loops. 



III. WANDERING LENGTH 

In one spatial dimension, the Pauli exclusion principle inhibits fermion permutations 
except in cases where the fermions wrap around the lattice boundary. For the remainder 
of our discussion, therefore, we consider systems with two or more dimensions. The first 
question we address is how far fermion worldlines can wander from start time t = to end 
time t = f3. We can put an upper bound on this wandering distance by considering the 
special case with no on-site potential and only nearest neighbor hopping. 

Let us consider motion in the x-direction. For each factor of S% in (^) a given fermion 
may remain at the same x value, move one lattice space to the left, or move one lattice space 
to the right. If h is the hopping parameter, then for large N the relative weights for these 
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possibilities are approximately 1 for remaining at the same x value, fthN' 1 for one move to 
the left, and fihN' 1 for one move to the right. In (|7|) we see that there are N factors of 
S^S®. Therefore for a typical worldline configuration at low filling fraction, /, we expect 
~ f3h hops to the left and ~ f3h hops to the right. For non- negligible / some of the hops 
are forbidden by the exclusion principle. Assuming random filling we expect ~ (1 — f)(3h 
hops to the left and ~ (1 — f)f3h hops to the right. 

The net displacement is equivalent to a random walk with 2(1 — f)(3h steps. The 
expected wandering length, A, is therefore given by 

A = yfiQ. ~ f)0h. (8) 

This result is somewhat surprising in that for typical simulation parameters (i.e., (3 not too 
large), we find A is no larger than a few lattice units. There is no contradiction between the 
existence of long distance correlations and the constraint of short distance wandering lengths. 
Long range signals are propagated by the net effect of many short range displacements. A 
simple analogy can be made with electrical conduction in a wire or sound propagation in 
a gas, which results from many short range displacements of individual electrons or gas 
molecules. 

In cases with on-site potentials, fermion hopping is dampened by differences in potential 
energy. Hence the estimate (J8|) serves as an upper bound for the general case. We have 
checked the upper bound numerically using simulation data generated by several different 
lattice Hamiltonians with and without on-site potentials. 

IV. PERMUTATION ZONE METHOD 

Let W be the logarithm of the partition function, 

W = log {Tr [exp(-/3#)]} . (9) 

Let us partition the spatial lattice, T, into zones Zi, Z 2 , Z^ such that the spatial dimen- 
sions of each zone are much greater than A. For notational convenience we define Z$ = 0. 
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For any R C V, let Wr be the logarithm of a restricted partition function that includes 
only worldline configurations where any worldline starting outside of R at t — returns to 
the same point at t — [3. In other words there are no permutations for worldlines starting 
outside of R. We note that Wr = W, and W$ is the logarithm of the restricted partition 
function with no worldline permutations at all. Since the zones are much larger than the 
length scale A, the worldline permutations in one zone has little or no effect on the worldline 
permutations in another zones. Therefore 

W ZoU .„ UZj - W Zo u...uz^ « W z . - W $ . (10) 
Using a telescoping series, we obtain 

W r -W d = £ (W Zo u...uz ] -W Zo u...uz^ 1 ) (11) 
j=\,...,k 

~ E (w Zj -w 9 ). 
j=i,...,k 

For translationally invariant systems tiled with congruent zones we find 

W = W r « W 9 + j§L(W Zl - W 9 ), (12) 

where |r|/|Zi| = k, the number of zones. For general zone shapes one can imagine 
partitioning the zones themselves into smaller congruent tiles. Therefore the result ( |P2"|) 
should hold for large arbitrarily shaped zones. For this case we take |T| to be the number of 
nearest neighbor bonds in the entire lattice and \Zi\ to be the number of nearest neighbor 
bonds in the zone. We will refer to \Zi\ as the zone size of Z\. This is just one choice 
for zone extrapolation. A more precise and complicated scheme could be devised which 
takes into account the circumscribed volume, number of included lattice points, and other 
geometric quantities. 

V. FREE FERMIONS 

As an example of the zone method, we compute the average energy (E) h^ 1 for a free 
fermion Hamiltonian with only hopping interactions on an 8 x 8 lattice. We consider values 
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(3h = 1.0, 1.5, and 2.0. The corresponding values for A are 1.0, 1.2, and 1.4 respectively. 
The Monte Carlo updates are performed using a single loop flip version of the loop algorithm 

In Fig. 3 we show data for rectangular zones with side dimensions 0x0, lx 1, 2 x 1, 
2x2, 3x2, 3x3, ...,6x6. We also show a least-squares fit (not including the smallest 



zones 0x0 and lxl) assuming linear dependence on zone size as predicted in (|12l). We 
find agreement at the 1% level or better when compared with the exact answers shown on 
the far right, which were computed using momentum-space decomposition. 

While the physics of the free hopping Hamiltonian is trivial, the computational problems 
are in fact maximally difficult. The severity of the sign problem can be measured in terms 
of the average sign, <Sign>, for contributions to the partition function. For (3h = 1.0, 
<Sign>~ 0.005; for (3h = 1.5, <Sign> ~ lO" 6 ; and for (3h = 2.0, <Sign>~ lO" 9 . Direct 
calculation using position-space Monte Carlo is impossible by several orders of magnitude 
for (3h > 1.5. 

VI. FERMIONIC 2D ISING MODEL 

While the free fermion example shows the linear dependence on zone size, we now study 
an example which better demonstrates the utility of the zone method. We will consider the 
Hamiltonian, 

H = [sijSi+xj + SijSij+x] (13) 

~h \ a i+X,j a i,j "I" a i,j a i+l,j a l,j+l a i,j a l,j a i,j+l 



h3 



where 



We will refer to this model as the fermionic 2D Ising model. When h = 0, our model 
reduces to the antiferromagnetic 2D Ising model for J > and the regular 2D Ising model 



for J < 0, with 2riij — 1 playing the role of Ising spin. We also note that the model 



for h > is equivalent to the h < model. On an even lattice, we can reverse the sign 
of h by multiplying —1 to all creation and annihilation operators on sites where i + j is 
odd. Unlike other quantum generalizations such as the transverse field Ising model (see for 
example JTT|] ) , we have reinterpreted the Ising spin as an occupation number and introduced 
nearest neighbor hopping of Fermi-Dirac particles. 

To our knowledge the fermionic Ising model has not previously been discussed in the 
literature. Therefore let us briefly discuss our interest in the model. The model is motivated 
by our studies of time-dependent background field fluctuations in Hamiltonian lattice gauge 
theories. In lattice gauge theory we encounter time-dependent Hamiltonians which in the 



Kogut-Susskind staggered formalism |L2| have the form 



H(t) = J2[H 1 ^(t) + H 2 , J (t)] + 



(15) 



where 



m,n 

//,,,(/) = E [-(^"-(*))*<] + i< - Km<hli + i 



(16) 
(17) 



The indices m and n are gauge group indices, and ti^^it) and hH£'™At) are background gauge 
fields. We will consider a simplified version of this system, where 

(18) 
(19) 



Hijj(t) — —hn j(t) ^aj +lj ajj + ajjaj+i^y . 
H2,i,j{t) = —h,2.i.j{t) {a\j +l aij + aj^aj 



and hi^jit) and ti2 1 i 1 j{t) are real valued functions. Let us consider what happens when 
hi,i,j(t) an d h 2t ij(t) are subject to Gaussian fluctuations. Let us define the Gaussian- 
integrated exponentials of the Hamiltonian, 



Ai ;i;j = j (II) . , exp [-dt ■ H 1;i;j (t)] exp -dt ■ ^j(h 1;i;j (t) - h) 



A 



%i,3 



dh 2 ,i,j(t) exp [-dt ■ H 2 ,ij(t)\ exp -dt ■ ^j(h 2 ^j(t) - h) 



(20) 
(21) 



We find 

A 1)itj oc exp 
A 2 ,i,j oc exp 

We note that 



dt • h jdij + aljtii+ij^ + cii • 2,7 (^J+ij + alja^ij^ 



(22) 
(23) 



(aS+ijOij + ajjfli+ij) = | (-Si+ijSij + 1) , (24) 

( a l,j+l a i,j + a ij a ij+l) — I ( — s i,j s i,j+l + 1) ) (25) 

and so the effect of such fluctuations is exactly modelled by an effective Hamiltonian of 
the form (|i~3D for J > 0. In the staggered formalism, fermion spin states are staggered at 
odd and even lattice sites, and the onset of antiferromagnetic order corresponds with the 
spontaneous breaking of chiral symmetry. 

The sign problem makes it difficult to study the fermionic Ising model at large volumes. 
There is no existing method that can handle this model with explicit fermionic degrees 
of freedom. For example near the transition temperature for h/J = 2.00 on a 20 x 20 
lattice we find <Sign>~ 10~ 4 . Given the significant computational difficulties we would 
like to see if the zone method could be effective for explicit fermion simulations near the 
critical temperature. One might expect the zone method to be useful since all of the 
physics of the h = Ising model is already contained in W$, the logarithm of the restricted 
partition function with no worldline permutations. For li ^ the contributions due to 
fermion permutations can be included in a controlled manner by considering successively 
larger permutation zones. 

Let us define the spontaneous staggered magnetization (S) as 

(S) = (-1)^' («y> = (-1)^' (2ny - 1) • (26) 

We will compute (S) for various temperatures and values h/J while adding a small external 
bias 



We compute (S) by linearly extrapolating results for different permutation zone sizes. For 
h/J< 1, the sign problem is significant but still manageable on a 20 x 20 lattice near the 
critical point. For h/J > 1, full simulations on a 20 x 20 lattice are not practical due to 
the sign. For h/ J = 1.50 we use permutation zones as large as 10 x 10, and for h/ J — 2.00 
we extrapolate using zones up to 6 x 6. In Fig. 4 we show the zone extrapolation for the 
staggered magnetization (S) at h/J = 1.50, T/J = 2.00, and m/J = 0.0256 on a 20 x 20 
lattice. In Fig. 5 we plot the zone-extrapolated results as a function oiT/J near the critical 
temperature for h/J — 1.50, m/J = 0.0256. Since the hopping interaction increases the 
disorder of the system, we expect the critical temperature to decrease as h increases. Noting 
that the h — Ising model has a critical temperature 

T c^U = i^y - 2.269, (27) 

we see that the data in Fig. 5 does in fact show a decrease in the critical temperature. We 
also find a critical exponent of f3 — 0.10(2), which is consistent with that of the Ising model, 
Asing = 0.125. 

In Fig. 6 we show the deviation of the critical temperature from T c Ising . In our plot we 
have included errors due to finite size effects, which we found by recalculating results using 
10 x 10 and 16 x 16 lattice sizes and several different values for m/J. We note that the 
deviation of the critical temperature appears to be quadratic in h/J, with possibly a small 
negative cubic contribution. It would be interesting to carry these simulations further to 
see how the curve behaves for larger values of h/J. About 2000 CPU hours on an IBM SP3 
were used to collect the data for this study. 

VII. SUMMARY 

We have presented a promising new approach to fermionic simulations at finite temper- 
atures. We have demonstrated that the exchange of identical fermions is short-ranged and 
has a maximum range of J 2(1 — f)f3h lattice sites, where (3 is the inverse temperature, h 
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is the hopping parameter, and / is the filling fraction. We have introduced the notion of 
permutation zones, special regions of the lattice where identical fermions may interchange 
and outside of which they may not. Using successively larger zones, one can extrapolate to 
obtain thermodynamic observables. We have demonstrated the method using two different 
examples: the average energy for a 2D free fermion and the spontaneous staggered spin 
magnetization for the 2D fermionic Ising model. 

The methods presented here are relatively easy to implement. Correlation functions 
can be calculated using worm algorithms ||13|| , which are generalizations of the loop algo- 
rithm. Although the full applicability of the zone method still needs to be determined, 
the method and its future generalizations may have some impact on several problems in 
strongly correlated solid-state systems and lattice gauge theory. In this paper we have 
focused on simulations with explicit fermion degrees of freedom. The extension of the zone 
method to pseudofermion algorithms is currently being studied and will presented in a future 
publication. 
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FIG. 1. Typical worldline configuration for the one-dimensional system. 



Arrow assignments 




FIG. 2. Upward/downward arrows are drawn at each occupied/unoccupied site. For each 
plaquette the number of inward-pointing arrows equals the number of outward-pointing arrows. 
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Zone extrapolations for <E>h~ l 
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FIG. 3. Average energy at j3h = 1.0, 1.5, 2.0 for free fermions on an 8 x 8 lattice. 
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FIG. 4. Zone extrapolation for the staggered magnetization (S) for the 2D fermionic Ising 
model at h/J = 1.50, T/J = 2.00, m/J = 0.0256. 
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<S> at h/J = 1 .50 for T <T. 
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FIG. 5. Staggered magnetization (S) for the 2D fermionic Ising model near the critical point 
for h/J = 1.50, m/J = 0.0256. 
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FIG. 6. Change in the critical temperature of the 2D fermionic Ising model as a function of h/J. 
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